# -*- coding: utf-8 -*-
'''
GNU GENERAL PUBLIC LICENSE
Version 2, June 1991
Copyright (C) 1989, 1991 Free Software Foundation, Inc.
51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
@author: Openagro
Creado en 08/11/2009
'''
import math

class Secador():
    """La clase del Secador es el objeto que se pretende modelar, el cual posee ciertos metodos
    y propiedades que se implementan a continuación"""
    def __init__(self,caudal,Tair,HRair,Tig,Hig,Hfg,msnm,Tamb,STS,ECG,int,intInv,HIp,TIp,ECGp):
        """Inicia la simulación llamando la función Thompson(), del modelo de Thompson.
        Ademas esta función es llamada para definir los datos como atributos de la clase"""
        self.caudal = caudal
        self.Tair = Tair
        self.HRair = HRair/100
        self.Tig = Tig
        self.Hig = Hig
        self.Hfg = Hfg
        self.msnm = msnm
        self.Tamb = Tamb
        self.STS = STS
        self.ECG = ECG
        self.int = int
        self.intInv = intInv
        self.HIp = HIp
        self.TIp = TIp
        self.ECGp = ECGp
        
        self.res_camSec = []
        self.res_camPreSec = []
        self.HIp_bs = self.chbh2chbs(self.HIp)
        self.HI_bs = self.chbh2chbs(self.Hig)
        self.res_camSec.insert(0, [self.HI_bs,self.HI_bs,self.HI_bs,self.HI_bs])
        self.res_camPreSec.insert(0, [self.HIp_bs,self.HIp_bs,self.HIp_bs,self.HIp_bs])
        
        
        #while True:
        for t in range(1,50):
            
            CHfg_bs1,CHfg_bs2,CHfg_bs3,CHfgprom,Tf3,Tf2,Tf1,CHfgprom_p = self.Thompson(self.caudal,self.Tair,self.HRair,
                          self.Tig,self.Hig,self.Hfg,self.msnm,self.Tamb,self.STS,self.ECG/3,t,self.intInv,
                          self.HIp,self.TIp,self.ECGp/3)
            
            self.res_camSec.insert(0,[CHfg_bs1,CHfg_bs2,CHfg_bs3,CHfgprom])
            
            CHfgprom_bh = self.chbs2chbh(CHfgprom)
            if CHfgprom_bh <= self.Hfg:
                break#Sale del While si se tiene una humedad promedio de las capas de grano menor  o
                #igual que la requerida 
            
            #Se calculan los parametros que se devolberan a la funcion Thompson
            #self.Tig = (Tf3 + Tf2 + Tf1)/3 #OJO No se enfria a lo sumo ley cero de la termodinamica
            self.Hig = self.chbs2chbh(CHfgprom)
            self.HIp = self.chbs2chbh(CHfgprom_p)
            #t = t+1
            #CHip,Tip,ECGp OJO no se hace presecado
        
    def Thompson(self,caudal,Tair,HRair,Tig,CHig,CHfg,msnm,Tamb,STS,ECG,t,intInv,CHip,Tip,ECGp):
        """Simulación del secado de cafe pergamino utilizando el modelo de Thompshon, secado en capa 
        Delgada"""
        
        #Condiciones sicrometricas del aire
        CHig_bs1 = self.chbh2chbs(CHig)
        Patm = self.PATM(msnm)
        Tair_abs = self.tmpC2K(Tair)
        Pvs = self.pvs(Tair_abs)
        PV = Pvs*HRair
        VE = (287*(Tair)/(Patm-PV)) #Volumen especifico del aire
        Vol = self.VolumenSecador(STS,ECG)#Volumen del silo
        RO = 365.884 + 2.7067*100*CHig_bs1 #Densidad aparente kg/m3
        ROS = RO*(1-CHig/100)#contenido de materia seca
        CAS = Vol*ROS#Capacidad de la camara de secado kg de cafe
        Qair = self.CaudalAir(caudal,STS)#Caudal en m3/hr
        R = VE*CAS/(3*Qair*1)#Relación grano aire
        
        
        #Inicio primera Capa
        cp1 = self.Cp(CHig,R)
        H1 = self.razH(PV,Patm)
        Te1 = self.bace(cp1,H1,Tair,Tig)
        Te1_abs = self.tmpC2K(Te1)
        Pvs = self.pvs(Te1_abs)
        HRair = self.calcHr(Patm, H1, Pvs)
        CHEQ = self.ROA(HRair,Te1)
        #DPV1 = Pvs1 - PV1
        DPV1 = Pvs*(1-HRair)/100
        dMdt1 = self.dm_dt(CHig_bs1,CHEQ,DPV1,t)
        CHfg_bs1 = CHig_bs1 + dMdt1
        Clv = self.clv(Tig,CHig_bs1/100)
        deltaH1 = (CHig_bs1/100 - CHfg_bs1/100)*R
        Hf1 = H1 + deltaH1
        Tf1 = self.bacs(H1,Te1,deltaH1,Clv,cp1,Hf1)
        #Fin primera Capa
        
        if Hf1 > 1:
            Hf1 = self.Condensacion(Hf1)
        
        #Inicio segunda Capa
        CHig_bs2 = CHfg_bs1
        CHig2 = self.chbs2chbh(CHig_bs2)
        Tig2 = Tf1
        Tair2 = Tf1
        Tair2_abs = self.tmpC2K(Tair2)
        Pvs2 = Pvs
        HRair2 = HRair
        PV2 = Pvs2*HRair2
        cp2 = self.Cp(CHig2,R)
        H2 =  Hf1
        Te2 = self.bace(cp2,H2,Tair2,Tig2)
        Te2_abs = self.tmpC2K(Te2)
        Pvs2 = self.pvs(Te2_abs)
        CHEQ2 = self.ROA(HRair2,Te2)
        DPV2 = Pvs2*(1-HRair2)/100
        dMdt2 = self.dm_dt(CHig_bs2,CHEQ2,DPV2,t)
        CHfg_bs2 = CHig_bs2 + dMdt2
        Clv = self.clv(Tig2,CHig_bs2/100)
        deltaH2 = (CHig_bs2/100 - CHfg_bs2/100)*R
        Hf2 = H2 + deltaH2
        Tf2 = self.bacs(H2,Te2,deltaH2,Clv,cp2,Hf2)
        #Fin segunda Capa
        
        if Hf2 > 1:
            Hf2 = self.Condensacion(Hf2)
        
        #Inicio tercera Capa
        CHig_bs3 = CHfg_bs2
        CHig3 = self.chbs2chbh(CHig_bs3)
        Tig3 = Tf2
        Tair3 = Tf2
        Tair3_abs = self.tmpC2K(Tair3)
        Pvs3 = Pvs2
        HRair3 = HRair2
        PV3 = Pvs3*HRair3
        cp3 = self.Cp(CHig3,R)
        H3 =  Hf2
        Te3 = self.bace(cp3,H3,Tair3,Tig3)
        Te3_abs = self.tmpC2K(Te3)
        Pvs3 = self.pvs(Te3_abs)
        CHEQ3 = self.ROA(HRair3,Te3)
        DPV3 = Pvs2*(1-HRair3)/100
        dMdt3 = self.dm_dt(CHig_bs3,CHEQ3,DPV3,t)
        CHfg_bs3 = CHig_bs3 + dMdt3
        Clv = self.clv(Tig3,CHig_bs3/100)
        deltaH3 = (CHig_bs3/100 - CHfg_bs3/100)*R
        Hf3 = H3 + deltaH3
        Tf3 = self.bacs(H3,Te3,deltaH3,Clv,cp3,Hf3)
        #Fin tercera Capa
        
        CHfgprom = self.ChfgProm(CHfg_bs1,CHfg_bs2,CHfg_bs3)
            
        if HRair3 > 1:
            HRair3 = self.Condensacion(HRair3)
       
        #Valida si la camara de secado ya llego al contenido de humedad deseado 
        
        
        #Una vez el aire cruza la masa de granos en la camara de secado se procede a 
        #simular el proceso de secado en la camara de presecado (Segundo piso del secador)
        
        CHfgprom_p = self.presecado(caudal,Tair,HRair,Tip,CHip,msnm,Tamb,STS,ECGp,t)
        
        #Validación del intervalo de inversión del flujo
        
        if t%7 == 0:
            if t == 1:
                pass
            else:
                CHfg_bs1,Tf1,HRair,CHfg_bs2,Tf2,HRair2,CHfg_bs3,Tf3,HRair3 = self.invFlujo(CHfg_bs1,Tf1,HRair,CHfg_bs2,Tf2,HRair2,CHfg_bs3,Tf3,HRair3)
        
        
        
        return CHfg_bs1,CHfg_bs2,CHfg_bs3,CHfgprom,Tf3,Tf2,Tf1,CHfgprom_p  #Fin del metodo Thompson()

    
    def ChfgProm(self,CHfg1,CHfg2,CHfg3):
        """calcula el contenido de humedad promedio CHfg"""
        CHfgprom = (CHfg1 + CHfg2 + CHfg3)/3
        return CHfgprom
         
    def presecado(self,caudal,Tair,HRair,Tig,CHig,msnm,Tamb,STS,ECG,t):
        """Metodo que implementa el ciclo de presecado"""
        #Inicia el calculo
        #Condiciones sicrometricas del aire

        CHig_bs1 = self.chbh2chbs(CHig)
        Patm = self.PATM(msnm)
        Tair_abs = self.tmpC2K(Tair)
        Pvs = self.pvs(Tair_abs)
        PV = Pvs*HRair
        VE = (287*(Tair)/(Patm-PV)) #Volumen especifico del aire
        Vol = self.VolumenSecador(STS,ECG)#Volumen del silo
        RO = 365.884 + 2.7067*100*CHig_bs1 #Densidad aparente kg/m3
        ROS = RO*(1-CHig/100)#contenido de materia seca
        CAS = Vol*ROS#Capacidad de la camara de secado kg de cafe
        Qair = self.CaudalAir(caudal,STS)#Caudal en m3/hr
        R = VE*CAS/(3*Qair*1)#Relación grano aire
        
        if t == 1:
            self.res_camPreSec.insert(0, [CHig_bs1,CHig_bs1,CHig_bs1,CHig_bs1])
            CHfgprom_p = CHig_bs1
        else:
            #Inicio primera Capa
            cp1 = self.Cp(CHig,R)
            H1 = self.razH(PV,Patm)
            Te1 = self.bace(cp1,H1,Tair,Tig)
            Te1_abs = self.tmpC2K(Te1)
            Pvs = self.pvs(Te1_abs)
            #HRair = self.calcHr(Patm, H1, Pvs)
            CHEQ = self.ROA(HRair,Te1)
            #DPV1 = Pvs1 - PV1
            DPV1 = Pvs*(1-HRair)/100
            dMdt1 = self.dm_dt(CHig_bs1,CHEQ,DPV1,t)
            CHfg_bs1 = CHig_bs1 + dMdt1
            Clv = self.clv(Tig,CHig_bs1/100)
            deltaH1 = (CHig_bs1/100 - CHfg_bs1/100)*R
            Hf1 = H1 + deltaH1
            Tf1 = self.bacs(H1,Te1,deltaH1,Clv,cp1,Hf1)
            #Fin primera Capa
            
            if HRair > 1:
                HRair = self.Condensacion(HRair)
            
            #Inicio segunda Capa
            CHig_bs2 = CHfg_bs1
            CHig2 = self.chbs2chbh(CHig_bs2)
            Tig2 = Tf1
            Tair2 = Tf1
            Tair2_abs = self.tmpC2K(Tair2)
            Pvs2 = Pvs
            HRair2 = HRair
            PV2 = Pvs2*HRair2
            cp2 = self.Cp(CHig2,R)
            H2 =  Hf1
            Te2 = self.bace(cp2,H2,Tair2,Tig2)
            Te2_abs = self.tmpC2K(Te2)
            Pvs2 = self.pvs(Te2_abs)
            CHEQ2 = self.ROA(HRair2,Te2)
            DPV2 = Pvs2*(1-HRair2)/100
            dMdt2 = self.dm_dt(CHig_bs2,CHEQ2,DPV2,t)
            CHfg_bs2 = CHig_bs2 + dMdt2
            Clv = self.clv(Tig2,CHig_bs2/100)
            deltaH2 = (CHig_bs2/100 - CHfg_bs2/100)*R
            Hf2 = H2 + deltaH2
            Tf2 = self.bacs(H2,Te2,deltaH2,Clv,cp2,Hf2)
            #Fin segunda Capa
            
            if HRair2 > 1:
                HRair2 = self.Condensacion(HRair2)
            
            #Inicio tercera Capa
            CHig_bs3 = CHfg_bs2
            CHig3 = self.chbs2chbh(CHig_bs3)
            Tig3 = Tf2
            Tair3 = Tf2
            Tair3_abs = self.tmpC2K(Tair3)
            Pvs3 = Pvs2
            HRair3 = HRair2
            PV3 = Pvs3*HRair3
            cp3 = self.Cp(CHig3,R)
            H3 =  Hf2
            Te3 = self.bace(cp3,H3,Tair3,Tig3)
            Te3_abs = self.tmpC2K(Te3)
            Pvs3 = self.pvs(Te3_abs)
            CHEQ3 = self.ROA(HRair3,Te3)
            DPV3 = Pvs2*(1-HRair3)/100
            dMdt3 = self.dm_dt(CHig_bs3,CHEQ3,DPV3,t)
            CHfg_bs3 = CHig_bs3 + dMdt3
            Clv = self.clv(Tig3,CHig_bs3/100)
            deltaH3 = (CHig_bs3/100 - CHfg_bs3/100)*R
            Hf3 = H3 + deltaH3
            Tf3 = self.bacs(H3,Te3,deltaH3,Clv,cp3,Hf3)
            #Fin tercera Capa
            
            if HRair3 > 1:
                HRair3 = self.Condensacion(HRair3)
            
            CHfgprom_p = self.ChfgProm(CHfg_bs1,CHfg_bs2,CHfg_bs3)
            dl = CHfgprom_p*0.02+t
            self.res_camPreSec.insert(0, [CHfg_bs1+dl,CHfg_bs2+dl,CHfg_bs3+dl,CHfgprom_p+dl])
        return CHfgprom_p

    
    def invFlujo(self,CHfg1,Tf1,Hr1,CHfg2,Tf2,Hr2,CHfg3,Tf3,Hr3):
        """Invierte las capas de grano en lacamara de secado, para recrear el fenomeno de inversión de flujo"""
        return CHfg3,Tf3,Hr3,CHfg2,Tf2,Hr2,CHfg1,Tf1,Hr1
        
    def Cp(self,CHg,R):
        """Calculo del calor especifico del grano en Kcal/(Kg de grano-C), Contenido de humedad %bs y la
        relación grano aire R kg mat seca de cafe/ kg de aire seco"""
        CHbs = self.chbh2chbs(CHg)
        cp = (1.3556 + 5.7859*CHbs)*0.239*R
        return cp
    
    def razH(self,PV,Patm):
        """Calculo de la razón de humedad Apartir de las ecuaciones sicrometricas"""
        H = 0.6219*PV/(Patm-PV)
        return H
    
    def ROA(self,Hr,Te):
        """Calcula el contenido de humedad en equilibrio para cafe pergamino 
        Utilizando el modelo desarrolado por Roa, apartir de Hr decimal, Teq en C, devuelve  CHEQ %bs"""
        P1 = 61.030848
        P2 = -108.37141
        P3 = 74.461059
        Q1 = -0.037049
        Q2 = 0.070114
        Q3 = -0.035177
        CHEQ = 0.01*(P1*Hr+P2*math.pow(Hr,2) +P3*math.pow(Hr,3))*math.exp((Q1*Hr+Q2*math.pow(Hr,2) +
                                                                         Q3*math.pow(Hr,3))*(Te))
        #CHEQ = CHEQ*100
        return CHEQ
    
    def clv(self,Tig,Chfg):
        """Calula el calor latente de vaporización, a partir de Temperatura inicial del grano en Celcius
         y el contenido de humedad del cafe en decimal base seca."""
        Clv = (25024 - 2.42958*Tig)*(1+1.44408*math.exp(-21.5011*Chfg))*0.239
        return Clv
    
    def bace(self,cp,H,Tair,Tig):
        """Balance de calor antes del secado"""
        Te = ((0.24 + 0.45*H)*Tair+cp*Tig)/(0.24+0.45*H+cp)
        return Te
    
    def bacs(self,Ho,Te,deltaH,Cvl,cp,Hf):
        """Balance de calor despues del secado"""
        Tf = ((0.24+0.45*Ho)*Te-deltaH*(587.9+Cvl-Te)+cp*Te)/(0.24+0.45*Hf+cp)
        return Tf
    
    def tmpC2K(self,T):
        """Conversión de temperatura de grados Celcius a grados Kelvin"""
        return T+273.16    
    
    def chbh2chbs(self,chbh):
        """Conversión entre contenido de humedad base humedad (chbh en %) y contenido de humedad 
        en base seca (chbs en %), retorna chbs en %"""
        chbs = chbh/(100-chbh)*100
        return chbs
    
    def chbs2chbh(self,chbs):
        """Conversión entre contenido de humedad base humedad (chbh en %) y contenido de humedad 
        en base seca (chbs en %), retorna chbh en %"""
        chbh = chbs/(100+chbs)*100
        return chbh
    
    def pvs(self,Tabs):
        """Calcula la presión de vapor de saturación apartir de Te (Temperatura de equilibrio: Balance antes del secado) en grados Kelvin"""
        VTRA = -7511.52/Tabs + 89.63121 + 0.02399897*Tabs - 0.000011654551 * math.pow(Tabs,2)
        VTRA = VTRA - 0.000000012810336*math.pow(Tabs,3) + 2.0998405e-11*math.pow(Tabs,4) - 12.150799*math.log(Tabs)
        Pvs = 1000*math.exp(VTRA)
        return Pvs
    
    def calcHr(self,Patm,H,Pvs):
        """Calcula la humedad relativa en fracción apartir de la presión atmosferica, 
        Razón de humedad y la presión de vapor de saturación."""
        Hr = Patm*H/(Pvs*(H+0.622))
        return Hr
        
    def Condensacion(self,Hf_r):
        """Rutina de condensación"""
        return Hf_r*0.9
    
    def dm_dt(self,chbs,CHEQ,DPV,t):
        """Ecuación de Secado en capa delgada"""
        m = 0.0143
        q = 0.87898
        n = 1.06439 
        w = q-1
        dMdt = -m*q*(chbs-CHEQ)*math.pow(DPV,n)*math.pow(t,w)
        return dMdt
    
    def PATM(self,msnm):
        """Calcula la presión atmosferica apartir de la altura sobre el nivel del mar msnm"""
        Patm = 9811.075*math.pow(10,1.014-msnm/19500)
        return Patm
    
    def CaudalAir(self,caudal,STS):
        """Calcula el caudal del aire"""
        Qair = caudal*STS*60
        return Qair
    
    def VolumenSecador(self,STS,ECG):
        """Calcula el volumen de grano en el secador"""
        Volumen = STS*ECG
        return Volumen
    
    def FlujoMasico(self,Qair,VE,STS):
        """Calculo de flujo masico"""
        fm = Qair/(VE*STS)
        return fm